iV-electron Slater determinants from 
non-unitary canonical transformations of fermion operators 



Carlos A. Jimenez-Hoyos, 1 R. Rodriguez-Guzman, 1 ' 2 and Gustavo E. Scuseria 1 ' 2 

1 Department of Chemistry, Rice University, Houston, TX 77005 
2 Department of Physics and Astronomy, Rice University, Houston, TX 77005 

(Dated: August 2, 2012) 

Abstract 

Mean-field methods such as Hartree-Fock (HF) or Hartree-Fock-Bogoliubov (HFB) constitute the build- 
ing blocks upon which more elaborate many-body theories are based on. The HF and HFB wavefunctions 
are built out of independent quasi-particles resulting from a unitary linear canonical transformation of the 
elementary fermion operators. Here, we discuss the possibility of allowing the HF transformation to be- 
come non-unitary. The properties of such HF vacua are discussed, as well as the evaluation of matrix 
elements among such states. We use a simple ansatz to demonstrate that a non-unitary transformation 
brings additional flexibility that can be exploited in variational approximations to many-fermion wavefunc- 
tions. The action of projection operators on non-unitary based HF states is also discussed and applied to 
the one-dimensional Hubbard model with periodic boundary conditions. 
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I. INTRODUCTION 



Mean-field methods such as Hartree-Fock (HF) or Hartree-Fock-Bogoliubov (HFB) have be- 
come paradigmatic in the description of many-fermion physics. These methods have found a wide 
range of applications in nuclear structure theory, condensed matter physics, and quantum chem- 
istry. This is not only because they constitute the simplest approximations to the exact many-body 
wavefunction, but also because more elaborate correlated approximations usually start from such 
independent quasi-particle vacua (HF or HFB). 

The HFB wavefunction developed to explain superconductivity relies on the so-called Bogoliubov- 
Valatin transformation, which defines quasi-particle operators as linear combinations of 

single-fermion creation and annihilation operators. These are then used to form a quasi-particle 
product state, the HFB wavefunction. Berezin studied the properties of general linear trans- 
formations of fermionic operators within a second-quantized framework. In this sense, one can 
consider the HF and HFB wavefunctions as being built out of single quasi-particle operators that 
result from a linear canonical transformation of the elementary fermion ones. 

A canonical transformation is understood in an algebraic framework as that which preserves 
the Dirac bracket of the phase-space variables in quantum mechanics (the position and momen- 
tum operators) [4]. In a second-quantized framework, this corresponds to a transformation that 
preserves the anti-commutation rules of the elementary fermion operators Ql. A linear canonical 
transformation does not need to be unitary, although Dirac [6] and Weyl 7|] showed that uni- 
tary transformations are canonical. Standard HF or HFB methods in several fields of many-body 
physics are usually carried out using a unitary canonical transformation. In this work, we study 
the possibility of constructing A^-particle Slater determinants resulting from non-unitary linear 
canonical transformations. The extension to HFB determinants will be discussed in a follow-up 
paper [8]. 

We note that non-unitary canonical transformations have been discussed in the literature before. 
They are discussed, for instance, by Blaizot and Ripka [5] in the general context of canonical 
transformations of second-quantized operators. They have been used by Balian and Berezin B]in 
the evaluation of matrix elements between two different Bogoliubov states. Zhang and Tang 101 ]. 
and later Ma and Zhang [111], have studied the properties of linear canonical transformations of 
fermion operators, including the non-unitary ones that we have just referred to. We also mention 
the work of Anderson |4|, where the properties of non-unitary canonical transformations have been 
discussed in a purely algebraic context, without reference to a Hilbert space. 



If a single Slater determinant is used as an ansatz for the many-fermion wavefunction, the full 
flexibility that a non-unitary canonical transformation affords is not evident because it does not add 
additional degrees of freedom to those existing in a unitary transformation. On the other hand, one 
can construct more general ansatze that use the flexibility of such a non-unitary transformation. 
We discuss here what may be the simplest, two-determinant ansatz that exploits all the degrees of 
freedom that define a non- unitary transformation for iV-particle Slater determinants. This idea has 
not been explored before in the literature. We here derive all expressions required for the evaluation 
of matrix elements between non-unitary based iV-particle Slater determinants. We also discuss the 
variational optimization of states based on a non-unitary HF-type canonical transformation. 

Our interest in non-unitary HF-type transformations originated from our recent work on pro- 
jected HF calculations for molecular systems 12], 131 and the two-dimensional Hubbard Hamil- 



tonian with periodic boundary conditions (PBC) 



14]. The idea of using a symmetry-projected 



HF state as an approximation to the many-body wavefunction was proposed by Lowdin 



15] as 



early as 1955. We, building on techniques developed and successfully applied in nuclear physics 
[5, 16-20], have shown that symmetry-projection out of the most general HF transformation yields 
a multi-reference type wavefunction which can account for a very significant part of the electron 
correlations. We have observed that, the more general the transformation we use (or the more 
symmetries that are broken), the better the resulting projected wavefunction is able to account for 
the correlation structure of the true Hamiltonian eigenvector. It is then natural to explore whether 
a non-unitary canonical transformation, which has more degrees of freedom than the unitary one 
commonly used, would yield additional flexibility for HF wavefunctions in general, and projected 
HF states in particular. This work describes our efforts along this line. We show that, indeed, 
using a non-unitary canonical transformation, one can build more flexible ansatze (based on N- 
particle Slater determinants) from which additional correlations can be accounted for in variational 
approximations. 

This paper is organized as follows. In section HH we discuss some general properties of linear 
canonical transformations of fermion operators. We proceed to show in section IIIII how to con- 
struct TV-particle Slater determinants based on such transformations. Section IIVI discusses our 
extension of Thouless' theorem for non-unitary Slater determinants. This is followed by section [VT 
where we use this theorem to derive the form of matrix elements between non-unitary TV-particle 
Slater determinants. In section IVH we introduce a two-determinant ansatz that displays the full 
flexibility of a non-unitary transformation. We show in section I VIII how such an ansatz can be 
used in projected HF approaches. This is followed by an illustrative application of the proposed 



wavefunction ansatze to the one-dimensional Hubbard Hamiltonian with PBC in section IVIIII 



II. CANONICAL TRANSFORMATIONS 



We start by introducing a set of fermion annihilation and creation operators c = {c^, ct}, which 



obey the standard anti-commutation relations 



0, 



0, 



Cfc,C 



(k\j) = S jk , 



where \k) ((k\) is a single-particle ket (bra) state. 

We now introduce a new set c 
one by the linear transformation 



We now introduce a new set of fermion operators f3 = {/3k, pt}, which is related to the original 




(1) 



where we have arranged the sets of fermion operators c and (3 into single columns. Here, U, V, Y, 
and X are arbitrary M x M matrices, where M is the dimension of the single-particle space. For 
compactness, we write the transformation defined by Eq. [T]as 



(3 = Tc. 



(2) 



It should be stressed that we have not enforced the relation ffl = {fly in Eq. Q] as this leads 
to a standard unitary transformation. One can show [sj] that the transformation is unitary if the 
matrix T satifies 



T* = aTa, 



(3) 



where the matrix a is given by 




(4) 



Here, Eq. [3] implies U = X and V = Y. 

We do insist, on the other hand, in making our transformation canonical, which implies pre- 
serving the appropriate anti-commutation relations, that is, 



0. 



0. 



It is not difficult to prove [5( that the transformation T is canonical if it obeys 

TaT T = a, 



(5) 



Using Eq. \5\ one can easily deduce the form of the inverse transformation 

fx V*\ 

T = I • (6) 




Equation [5] also provides the conditions that the matrices U, V, X, and Y must satisfy for T to 
define a canonical transformation. Those are given by 

tri" x + yty = 1, (7a) 

X T U* + Y T V* = 1, (7b) 
[/t y* + yt J/* = 0j (7c) 

Y T X + X T Y = 0. (7d) 

Note that the matrices U> V* and Y T X are anti-symmetric. 

The matrices T form a group (the fermion group described by Ma and Zhang [llj]) isomorphic 
to the group of orthogonal matrices of dimension 2M [0(2M, C)] [5]. On the other hand, the set 
of matrices T for which the transformation is unitary form a group isomorphic to the group of real 
orthogonal matrices of dimension 2M [0(2M)j. There are twice as many degrees of freedom in 
choosing a general non-unitary transformation than in a unitary one. 

We close this section by noting that the transformation defined in Eq. Q] is more naturally 
understood as a linear transformation if one introduces an operator S such that 

(3 = ScS- 1 = Tc. (8) 



The form of the operator S has been discussed by Blaizot and Ripka [5], Zhang and Tang [10], and 



Ma and Zhang 



111- 



B. 



III. iV-ELECTRON SLATER DETERMINANTS 

In this section, we discuss the construction of iV-particle Slater determinants using quasi-particle 
operators resulting from canonical transformations of the elementary fermion ones. This is dis- 
cussed in detail by Navon [21[, as well as in several textbooks in many-body physics. 



In standard (i.e. unitary) HF theory, an iV-electron Slater determinant is constructed out of a 
set ./V hole creation ({fejj}) and M — N particle annihilation ({b p }) operators, each of them resulting 



from a linear combination of the elementary operators {c k , cl}: 

3 

b p = Dj p Cj . 



(9a) 
(9b) 



Using standard notation, the first ./V columns in D (which we write as Dh) represent the hole 
states, while the last M — N columns (which we write as D p ) represent the particle states. 

The transformation from the elementary operators to the set of HF operators constructed above 
can be written as 







b p 




b h 




\bi) 





'NxM 



\ 







(M-N)xM 







NxM 



xM 



D 



(10) 



/ 



where we have implicitly assumed the transformation to be unitary. 

The above transformation is canonical if the HF operators satisfy the (non-trivial) anti- 
commutation relations 



b h , b{, 



b p , &J, 



5 pip, 



b p ,b[ 



0. 



These conditions restrict the form of the matrix D according to 



jk 



6 P' b l'\ + = Yl D iv D l P > Sjk = [D ] D) nln = Vp> 

jk 



P P 



b p ,b[ 



jk 



(11a) 
(lib) 
(11c) 



The first equation implies orthonormality of the hole states, the second one orthonormality of the 
particle states, and the last one corresponds to orthogonality between hole and particle states. All 
these conditions are summarized in the requirement D = 1. 

One could allow the HF transformation described previously to become non-unitary by intro- 
ducing, in addition to the operators described by Eq. [91 another set of hole and particle operators, 



{bh,V p }, given by 



j 



(12a) 
(12b) 



A non-unitary transformation can then be built as 

' OnxM D 



h 
\P P J 



t \ 



Dj 



(M-N)xM 

Onxm 

\®(M-N)xM Dp J 



It is a canonical transformation if the (non-trivial) anti-commutation relations 



(13) 



b p ,b\ 



0. 



bh,bp 



0, 



bh,bl, 



$h'h, 



bp, b\, 



8 'pip, 



are satisfied. These conditions restrict the form of the matrices D and D according to 



v b l] + = E ^> s,k = D ) hp = °> 

jk 



h, bp] = Y, D 3h Dip S jk = (5t D) = 0, 

J + '—f c \ / ph 



bh,bl 



Y J D jh Dl h ,6 jk =[rtD) hih = 8 h 



bp, bl 



Y,D jp D* w 5 jk =(D^D)^ = 8, 



j Pp. 



(14a) 
(14b) 
(14c) 
(14d) 



The first two equations imply orthogonality of the hole and particle states in D and D. The last 
two equations imply a bi-orthonormality between the hole and particle orbitals in D and D. Note 
that the last two conditions are satisfied by choosing D< = D^ 1 , but the orthogonality among hole 
and particle states has to be separately imposed. 

Let us remark that, if the HF operators {b^, b p , b~h, bp} define a canonical transformation, the 
inverse transformation is given by (see Eq. [6]) 







u 


l-( 



OmxJV D* 



Dfr ®Mx{M-N) 



Dh Omx{m-n) OmxN Dp 



h 

\blj 



(15) 



The bi-orthonormal Slater determinants |<3?) and |$) are produced when the set of operators 
{6^,6^} act on the bare fermion vacuum |— }, i.e., 



i$> = IKi->' ( 16 ) 

h 

i*>=IKi->- ( 17 ) 

h 

They satisfy the bi-orthonormality condition (^l^) = 1. 

One can easily show that \&) and |<3?) act as vacua to a certain set of hole or particle states: 

b{\<S>) =0 V b[, b p \<S>) = V b p , 
bl\$)=0 V b[, 6 P |¥) = V b p . 

IV. THOULESS' THEOREM FOR TV-ELECTRON SLATER DETERMINANTS 

In standard (i.e. unitary) HF, there is a theorem due to Thouless [221 ] which reads: 

Theorem. Given a Slater determinant |<3?o) which is a vacuum to the operators 
{bh,b p }, any iV-particle Slater determinant |$i) which is not orthogonal to l^o) can 
be written in the form 

=A/"exp (^Zphblb^j |$ ), (18) 

where M = (^ol^i) is a normalization constant and the coefficients Z v h are uniquely 
determined. Conversely, any wavefunction of the form of Eq. [T51 where |$o) is a Slater 
determinant, is also an iV-particle Slater determinant. 

For Slater determinants built out of operators resulting from a non-unitary linear canonical 
transformation, the equivalent theorem reads 

Theorem. Given a Slater determinant l^o) which is a vacuum to the operators 
{b^bp}, any iV-particle Slater determinant which is not orthogonal to j^o) can 
be written in the form 

|* x ) =A/"exp {Y,Zp h blb^\ |$o), (19) 

where M = (^ol^i) is a normalization constant and the coefficients Z p h are uniquely 
determined. 

For a proof of this last theorem we refer the reader to Appendix [A] of the present work. 



V. MATRIX ELEMENTS BETWEEN iV-ELECTRON SLATER DETERMINANTS 



In this section we obtain the expressions required for the evaluation of matrix elements be- 
tween arbitrary Slater determinants built out of operators resulting from a non-unitary canonical 
transformation. 



A. Norm overlaps 



The overlap between two iV-particle Slater determinants of the form |$ a ) = flfc !! - ) can De 



obtained by application of Wick's theorem 



5] 



on the bare fermion vacuum. That is, 



{§p\$ a ) = H/3;v-/Wi---4r|-) =det5, (20) 

where Sij = ftaj = (/3j|ay). Here, we have used the fact that the contractions fli (3j and a\ 
vanish for HF-type operators. 

The overlaps among TV-particle Slater determinants become 

($o|*i) = detjv D 0T D 1 *, (21a) 

($ |li) = detjy D 0T D 1 *, (21b) 

($ |$i) = det N D 0T D 1 *, (21c) 

($ |$i) = detjy £> 0T (21d) 

where we have used detjv to denote that the determinant is over the N x N set of occupied or- 
bitals. Observe that (^ol'&o) = (^ol^o) = 1> which corresponds to the bi-orthonormality condition 
previously described. 



B. Operator matrix elements 

In deriving the expressions for operator matrix elements, we follow Ring and Shuck [16]. Our 
aim in this subsection is to evaluate matrix elements of the form 

($o\4 1 ---4 p Ck 1 ---c kp \$i). (22) 

The form above is chosen for convenience, but other matrix elements can be derived in the same 
way described below. 



We shall use Thouless' theorem to write the state |3>i) as 

|* 1 } = exp(i)|$ ){¥o|$i}, (23) 
Z = Y,Z P h%b h . (24) 

ph 

Here, {b h ,b p ,bh,b P } are defined such that 

b{\$ )=0 V b[, b P \$ } = V b p , 
(¥ |6 ft = V b h , (*o|3 = V 4- 

On the other hand, we write the state (<J?o| as 

(¥ | = (¥ |exp(-i), (25) 

where use has been made of the vacuum properties just described. 

It then follows that we can evaluate the general matrix element from Eq. [22] as 

(®o\cl ■■■c\ p c kl --- c kp \<$>i) = ($ |$l)(!o| exp(-i) cl ■ ■ ■ c\ p c kl ■ ■ ■ c kp exp(i)|$ ) 

= (¥ |$i)(¥ |4 • • • d lp d kl --- 4J$ ), (26) 

where we have introduced the operators 

d\ = exp{—Z) c\ exp(Z), (27a) 
d k = exp(-Z) c k exp(Z). (27b) 

We now express the operators {d~i, dk} in terms of {b k , b p , bh,b P }. This is accomplished by using 
Eq. [15]to write {cj,cj} in terms of {b h , b p , bh, bp}- It follows that 

d\ = exp(— Z) c\ exp(iJ) = c\ — \z, c\ 

= E A , ^i + E ( D l - E 3* A ,) % , (28) 

h p V h J 

d k = exp(-i^) c k exp(i^) = c k - \z, c fc J 

= E ( D lh + E Z Ph D k*pj h + E D t V (29) 

h \ P / P 

Because {d~i,dk} are given as linear combinations of {b k , b p , bh, bj,}, Wick's theorem [5j can be 



used to calculate the corresponding matrix elements. The non- vanishing contractions among the 

10 



operators {b^, b p , bh, b P } are given by 

K b h' = S hh', (30a) 
bpV pl = 5pp>. (30b) 

It follows that the non- vanishing contractions among the operators are of the form 

did k = } j D? h ( Dl* h , + Z p/l / Z)£* j 
Wi' V p / 

= £ A°. Dll + £ D,° h ^ 5g, (31) 

h ph 



S fe = 53 A ; IK* - E z ^ D kh) 

pp' \ h / 



Sppi 



= ^ A? D° kp - ^ Ap ^ D° kh . (32) 
p pfe 

The application of Wick's theorem to the operator matrix elements of the form of Eq. [22] leads 
us to conclude that all such matrix elements can be evaluated in terms of the transition density 
matrix p 01 , given by 

p oi = = ($ |exp(-i) C t Cfc exp(i)|d> ) 

\*o|*i/ 

= ^ Df h D° k * h + 53 Dl Z ph Dtp, (33) 

h ph 

where 

£ »=( j5o,j,, ) M - < 35 » 

Here, we have used Eqs. I A2I and I A5I from Appendix|A]to write the forms of the matrices Z and C 



C. Evaluation of the energy of a single Slater determinant 

As an example of the application of the above equations, let us now consider the evaluation 
of the energy of a determinant |<&). Given a two-body Hamiltonian in the usual second-quantized 
form l^l 

h = 53^1^) 4 °k + 1 ^2{iM*>i) 4 c) a ok, (36) 

ik ijkl 
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where (i\h\k) and (ij\v\kl) are one- and anti-symmetrized two-particle integrals, respectively, the 
energy can be evaluated as 



where 



and 



E 



(*l*> 



2 

Tr^p+irpV (37) 



($|ctc fc |$) 

= E a* ^ + E ^ ( 38 ) 

r ifc = E(^>l fc 0wi, (39) 



Z ph = ^ (D J D* )(C*- l ) h , h , (40) 
h> p 



C h , h =[tfD) . (41) 
\ J h'h 

It is important to realize that the energy expression (Eq. |37|) has the same form as in standard 
(i.e. unitary) HF. The difference lies in the form of the density matrix p (Eq. |38|) . which comes 
about from the fact that the anti-commutation relations satisfied by the HF operators are different. 



VI. VARIATIONAL ANSATZ WITH SLATER DETERMINANTS FROM NON-UNITARY 
TRANSFORMATIONS 

In this section, we use a simple, two-determinant ansatz that uses the full flexibility of the 
non-unitary HF-like transformation of Eq. [13] as part of a variational strategy. 

Before introducing such ansatz, we note that using a single Slater determinant |$) as a trial 
wavefunction, whether resulting from a unitary or a non-unitary canonical transformation, would 
lead to the same variational energy. An JV-particle Slater determinant resulting from a non-unitary 
canonical transformation is equivalent to an un-normalized Slater determinant in the usual (i.e. 
unitary) sense. The variational optimization of the energy (taken as the Hamiltonian overlap over 
the norm overlap) would lead to the same result regardless of the underlying normalization of the 
determinant. 
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The two-determinant ansatz that we use is given by 

|*) = ci|$) + c 2 |$), 

= ci|$i) + c 2 |$ 2 ), (42) 

where c\ and c 2 are coefficients to be determined variationally. We have made the identification 
\<&i) = |<3?) and |<I> 2 } = 1$} to simplify our notation below. Observe that for a standard {i.e. 
unitary) HF transformation, = |$»), which in turn implies l^} = |<3?). 

One could argue that the ansatz of Eq. [32] has the same variational flexibility as that in which 
|<J>i) and |<I> 2 ) are two non-orthogonal Slater determinants resulting, each of them, from a standard 

we use explicitly results 
from a single linear canonical transformation of the elementary fermion operators. 

The Hamiltonian expectation value associated with the state |^) is given by 

2 



£ c* Q cp^ a \H\^ p ) 



E = ^ . (43) 

a,/3=l 

We rewrite the energy above in the form 

y a p = — • (45) 

C*'C^($ a '|^') 

a',/3'=l 

The matrix elements appearing in Eqs. EH and 05] can be evaluated in a straight-forward way. 
The overlap kernels in Eq. 05] are computed as 

($|$) = det N D T D*, (46a) 

(3F|S) = detArZ) 1 "!)* = 1, (46b) 

= detivD T i5* = 1, (46c) 

($|$) = det7v£> T i5*. (46d) 

The Hamiltonian kernels are evaluated in terms of transition density matrices as 

r? = £<<j|«|w>pjf. (48) 
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The transition density matrices are in turn given by 



\ 


&.{<&) 








/at « 1 (f>\ 




(*[*> 


<* 


4^1$) 


t 











Here, 



Pfci = 7^ = L D * D *kh + L ^ ( 49a ) 

h ph 

pJN^^-^A^, (49c) 



^ = \^ 1 = Ah D* kh + ^ D ih Z ph D* kp . (49d) 

h ph 



Z ph = Y J {n T D*) phi {C*- 1 ) h , h , (50a) 



h> 



^=EK 5 1„,( r_1 )'.".' (50b) 

£ h , A = (£)t £,) , (50c) 

V / h'h 

C h , h = D) . (50d) 
A. Variational optimization of 1$) 

Let us now consider the variational optimization of the wavefunction ansatz introduced in Eq. 
The variational parameters are the coefficients c\ and C2 and the orbital coefficients (that is, 
the matrices D and D) defining the states |$) and |<J>). The variation has to be carried out subject 
to the constraint that (^|^) = 1, which is equivalent to saying that |<£) and |$) are defined by a 
canonical transformation of the form of Eq. [T3l 

The variation with respect to the coefficients c\ and C2 yields the generalized eigenvalue problem 

(H-M)c = 0, (51) 

with the constraint 

C tNc = l, (52) 

which ensures the orthonormality of the solution. Here, c represents the column of coefficients 
{01,02} , while H and N are, respectively, Hamiltonian and overlap matrices given by 

Hap = {* a \H\$p), (53) 
N a0 = ($a\*p). (54) 
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It should be stressed that at this level we only keep the lowest-energy solution to the generalized 
eigenvalue problem, in a similar way as in projected-HF methods involving an eigenvalue problem 



Let us now consider the variation in the energy with respect to the underlying non-unitary HF 
transformation. We have followed the work of Egido and coworkers [23] for this purpose. Let us 
assume that we are provided a guess for |$) and |$), characterized by the set of HF operators 
{b h ,b p ,i>h,bp}. We can now parametrize the energy functional around {|$), l^)} by allowing for 
independent Thouless' rotations of both states, characterized by the matrices Z and Z. That is, 
we let 

|$) -> exp Z ph P p b h |$), (55a) 
\ P h J 

l ¥ > -> ex P f zZ 2 vh bt b h J | $) . (55b) 



We define the local gradient {G, G} around Z = and Z = as 



d 

G ph = ~ q-^T E l Z > Z \ 

d 



(56a) 
(56b) 



Here, Z p h and Z* h are treated as independent variables, and the same is true for Z p h and Z* h . The 
total derivative of the energy then becomes 

dE = ~ 2 [°ph dZ * P h + G ph dZ* ph + c.c] . (57) 

ph 

Explicit differentiation of the parametrized energy functional leads to the following expressions 
for the local gradient: 

($\P h b p (h-e) |$> (^\b[bp (h-e) |¥) 

G ^ = — wr^ — m — m (58a) 

(¥|6+ bp (6 -e)\$) @\bt bp (6 - e) |$) 



h P \ ^ \^\ u h U P . 

($1$) ($1$) 

where E is the energy corresponding to the state |^) from Eq. 

15 



The overlap-like matrix elements appearing in Eq. [58] can be evaluated as 



^ ] -^m/i -^np Prim' (59a) 
mn 

0, (59b) 

0, (59c) 



7=j=7 - 2^ rnh U np P nm - \p\)a) 

The Hamiltonian-like matrix elements appearing in Eq. [58] can be evaluated as 

' 1 ' mn ' 1 ' 

+ E E + r ^) ^™ ( 5 ™ ~ ^) ' ( 60a ) 

mn ik 

J2 D *khD ip {h ik + T%), (60b) 

ik 

J2 D thD ip {hk + Tli), (60c) 



(*[*> 



($|6[6 p ff|$; 



tfc 



/$|$\ ~ 2^ U rnh^n P P nm ,™ 

> 1 ' mn < 1 ' 

+ E E ^ ^ + P i) Pf™ (Sni ~ (%) ■ (60d) 

mn jfe 

B. Restoration of the bi-orthonormality condition 

Let us assume that, during the optimization process, we started with the states |<J?) and |3>) 
and produced the new states [$') and |$') by using the Thouless' transformations 

|$') = A/"exp (j2z ph blbl\ |$), (61a) 

II 7 ) = N exp ^ Z p/i 6j | $) . (61b) 

Here, the matrices Z and Z can be chosen as, for instance, 

Z p h = vGph, (62a) 
Z ph = 7]G ph , (62b) 
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with 7] > being some parameter. We denote with {d k , d p , dh,dp} the set of HF operators produced 
by such transformations (see Eqs. IA6al and IA6bj) 



A = b{ + Y. z v*~ h l 

p 

dp = bp — Z ph b h 

h 

dh = b h + Yl ^ph bp 

p 

S p = bl-^z; h b{ 



(63a) 
(63b) 
(63c) 
(63d) 



where the operators {b' h , b p , b^, b p } describing the states |$) and |$) are assumed to satisfy all the 
appropriate anti-commutation relations. 

We show in Appendix |A] that the operators {d h ,d p } annihilate the vacuum |$'}. Similarly, 
the operators {d\,d p } annihilate the vacuum |<fr'). The operators {d k , d p , dh, dp} do not, however, 
satisfy the anti-commutation relations given by Eq. Q3J In fact, they satisfy 



d rift 
u p , u h 



dh,dp 

4,4' 
dp, 4/ 



I + Z T Z* 



I + Z*Z T 



h'h 



VP 



(64a) 
(64b) 
(64c) 
(64d) 



We can restore the desired anti-commutation relations by performing the transformations 



4 - y^. L hh> d h >i 

h' 

4 = ^l* h -}d hl , 
h' 

4 = E M ;/4', 
P > 

p> 

in terms of the lower triangular matrices L, L, M, and M 



(65a) 
(65b) 
(65c) 
(65d) 



23]. 



The anti-commutation relations among {dl, d p , dh, dp} become 



4,41 =E z v lL ^ (i+z T z* 



fJ,V 



41 = E m p^ u p^ ( 7 + 2 * zT 



[MP 



UfJ, 



•'P'p, 



(66a) 
(66b) 
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which yield the following equations for determining L, L, M, and M: 



I + Z T Z* =LL\ 



I + Z* z x 



MM ] . 



(67a) 
(67b) 



Hence, given the matrices Z and Z, standard LU decompositions (Eqs. I67al and I67bp can be 
performed to obtain the matrices L, L, M, and M. This is similar to the unitary case, where the 
only two matrices required (L and M) can be obtained by Cholesky decompositions fl4. 

We remark that if Z = (or Z = 0), then the operators {d/ h , d p , d~h, dp} do obey all the required 
anti-commutation relations. In other words, one has to restore the bi-orthonormality condition 
only if both |<I>) and |$) are rotated. 



C. Global gradient 



In order to use gradient-based optimization methods such as the conjugate gradient or quasi- 



Newton methods (see Refs. 



14, 



20 



23, and 



24 ). one must be able to compute a global gradient. 



That is, we should be able to compute the gradient of the energy at with respect to variations 
in Z and Z defined in terms of th e op erators {b^, b p , b° h , Vp} corresponding to the reference state 
l^o}- Here, we follow Egido et al. 23j in deriving the form of the global gradient. 
Consider the energy of the state |\E'i). It is given by 



c* a cp^l\H\^}) 

a,p=l 
2 

a,P=l 



(68) 



Provided that l^ 1 } and l^ 1 ) are non-orthogonal to (<3?°| and (3>°|, respectively, we can write 



ph 



id> 1 )=Arex P l^2z ph b o p n° h U }, 



ph 



(69a) 



(69b) 
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where J\f = ( < I >0 |<I> 1 } and N = ( < I )0 | < J )X ) are normalization constants. Here, 



£>h'h 
£h'h 



h'h 



h'h 



(70a) 

(70b) 

(70c) 
(70d) 



where we have used Eqs. IA2l andl A5l to write Z and Z in terms of the matrices of orbital coefficients 
D°, D°, D 1 , and D 1 . 

A variation in Z and Z leads to a change in energy given by 



8E , d£ 

■ oZ nh + - - „ oz„ 



+ c.c. 



pll 

where we have introduced the global gradients Q and Q given by 



(71) 



Gph = —yu 



2/12 



(^Ij* 1 ) 



£ph = -2/21 , _ : , 2/22 



$ $1 



(72a) 



(72b) 



In order to evaluate the matrix elements appearing in the global gradient (Eq. I72|) . we need to 
relate the operators {b^, bp, W h , Vp} to the operators {b^ , bp, b^, bp* }. Combining the results of the 
previous subsection with Eqs. IA6al and IA6bl we arrive at 



# = E ^ © = E ^ U + E ^ ^1 . 

h' h' \ P / 

bl = E ^ % = E ^ (# - E ^ ^ 
bl = E z w ^ = E ^ ( + E %v ■ 

h' h' \ P / 

9 = E ^ = E (V - E ^ #1 

p' p' \ h / 



(73a) 
(73b) 
(73c) 
(73d) 



where the matrices L, L, M, and M are here determined by the solution to Eqs. I67al and I67bl 
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Because the transformation denned by Eqs. !73aH73dl is canonical (we have explicitly ensured 
that anti-commutation rules are preserved), we can invert the transformation using Eq. as a 
reference. We arrive at 



h' pp' 

^ = E^J' + E^' L W^> 

p' hh' 

= E L h^h bh' - E ^P'h ^PP 1 & P' 
hi pp' 



(74a) 
(74b) 
(74c) 
(74d) 



hh' 



We now use Eqs. !74aH74dl to write the global gradient {Q and Q) matrix elements in terms of 
the local gradient (G and G) as 



Gph = E L h'h M p'l G P'h' 



M T-lQ L *-l 



M T ~ 1 GL*- 1 



ph 



ph 



(75a) 
(75b) 



Gph ~ E L h'h M p'l G V'< 
p'h' 

We close this subsection by noting that one has reached a solution to the variational equations 
when the local gradient (and, consequently, the global gradient) vanishes, i.e., 

a 



az. 



E = Q, 



ph 







az 



£7 = 0. 



(76a) 
(76b) 



ph 



VII. VARIATIONAL ANSATZ WITH PROJECTION OPERATORS 



We now turn our attention to states resulting from the action of symmetry-restoring projection 
operators on symmetry-broken determinants. We start by providing a brief description of the form 

s.y, 
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18. 



of the projection operators used. More details can be found in Refs. 

Consider a symmetry group G, with elements {<?}, that commutes with the Hamiltonian. The 
group can be continuous or discrete, but we shall assume for simplicity that it is Abelian. A Slater 
determinant is symmetry broken if 



(77) 



that is, if the determinant is not invariant upon action by the elements {<?}. The set of all {^l^}} 
is called the Goldstone manifold. The norm and the matrix elements of commuting observables are 
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25] 



the same within the Goldstone manifold up to an arbitrary phase factor |5|. It is well known 
that the symmetry can be restored by diagonalization of the Hamiltonian among the Goldstone 
manifold. 

A projection operator can, in general, be written as 

pi = I f d6w J (6)R e , (78) 

where L is the volume of integration, Rq is an element of the symmetry group in consideration, 
the index j labels the eigenvalue restored by means of the projection, and the coefficients w J (9) 
correspond to the matrix elements of the operator Rq among the irreducible representations of the 
group. Evidently, for discrete groups the integration above is replaced by a discrete sum. We shall 
drop the label j henceforth for simplicity of notation. 

As an example of the projection operators discussed above, S z projection on a broken-symmetry 
determinant can be accomplished by 

P m = — [ dO exp \i0 (§ z -m)], (79) 



4vr „ 

where an eigenfunction of S z with eigenvalue m is recovered upon the action of the projection 
operator above. 

We work with cases where Rq are single-particle rotation operators that act on the HF ones 
according to 

b\(9) = Re blR, 1 = £ D* k Rq c )kf = ^ R^O) D* k cj, (80) 

i ij 

where Rij{6) = [i\Rg\j) is the matrix representation of Rq in the single-particle basis. 

We can now use the variational ansatz introduced in Eq. 02] and put a projection operator in 
front of it. The proposed wavefunction becomes 

P\V) = Jd9w{9) dRe\^) + c 2 Rq\$) . (81) 

The Hamiltonian expectation value of a wavefunction of the form of Eq. [81] can be written as 

= (*\piHP\*) = imPm 

ya/3{V) - — g , (83 ) 

d9w(9) cZ,cp($ Q ,\Rg\$ fi ,) 

a',/3'=l 
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where we have made the identifications |3>i) = |<5) and \& 2 ) = \&)- The expressions for the matrix 
elements appearing in Eqs. 1821 and 1851 are given in Appendix [Bj 

A. Optimization of the projected ansatz 

Our task is now to minimize the energy of our ansatz for the projected state (Eq. I82|) with 
respect to variations in the reference determinants |$) and |<I>). We will closely follow the derivation 
we presented before (section IVI Ap for the optimization of the unprojected state. 

The variation with respect to the coefficients c\ and c 2 yields a generalized eigenvalue problem 
similar to the one of Eqs. [5TJ and [52j In this case, H and N are 2x2 matrices given by 

H af) = J de w {6) (<f> a \HRe\$p), (84) 

N aP = J Mw{p){* a \Re\*f,). (85) 

Once again, only the lowest-energy solution is used in the variational optimization. 

The parametrization of the energy functional with respect to the determinants |$) and |<J>) is 
done in the same way as it was done for the unprojected case [l^Q]- That is, we parametrize the 
energy functional in terms of the Thouless' rotation matrices Z and Z acting upon |$) and |3>), 
respectively. 

The resulting local gradient is derived by using the definitions in Eqs. I56al and I56b[ We arrive 
at the expressions 



G ph = J d6w(9) < 



<$ fit b p (H - E) Rg\$) (* 6+ b p (H - E) R e m 
-yu{6) ^ yn(e) V s J V, (86a) 

mb[b p (h-e) R e m mb[b p (h-e\ R e m 



G ph = I d9w(6) { -y 21 {6) h J y 22 (9) * P _L ' — ) . (86b) 

Here, E is the energy corresponding to the state \^) from Eq. [8T1 The explicit expressions for the 
matrix elements appearing in Eq. [86] are given as part of Appendix [Bj We finally note that the 
relationship between the local gradient and the global gradient is the same as in the unprojected 
case (see Eq. 1751) . 

VIII. APPLICATION TO THE ONE-DIMENSIONAL HUBBARD HAMILTONIAN 



In this section we present the application of the ansatze discussed previously to the one- 
dimensional Hubbard Hamiltonian 26j with PBC. This describes a set of electrons in a lattice 
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according to 



Here, cj CT creates an electron on site j of the lattice with a = ^-projection of spin. The first 

term in the Hamiltonian accounts for a negative (t > 0) kinetic energy that the electrons gain when 
they hop from one site to a neighbor. The second term accounts for the (U > 0) repulsion that 
opposite-spin electrons feel when they are in the same site. The lattice used for this Hamiltonian is 
a finite one with N s sites. Periodic boundary conditions are assumed, which make the site N s + k 
equivalent to the site k. 

The ID Hubbard Hamiltonian has been extensively studied, and our purpose here is merely to 
test the flexibility that TV-particle Slater determinants constructed in terms of non-unitary canonical 
transformations bring. With this in mind, ours should be regarded as a proof of feasibility for 
calculations in finite many-fermion systems based on non-unitary HF transformations. We point 
the interested reader to the comprehensive book on the ID Hubbard Hamiltonian by Essler et al. 



261 ] . Recent work on the ID Hubbard Hamiltonian with projected HF approximations has been 



281 ] devised a set 



done by Schmid et al. [2j] and Tomita [271 ]. We also note that Lieb and Wu 
of equations from which the exact eigenvalues of the ID Hubbard Hamiltonian of Eq. [87] can be 
obtained. 

We have applied the methods described in the preceeding sections to the lD-Hubbard Hamil- 
tonian. Our calculations have been performed with an in-house code using the conjugate- gradient 
method described here and in Ref. [23j for the variational optimization of HF-based states (see also 
Refs. [2J and Il4). We have selected U = At as a representative on-site repulsion, corresponding 
to a strongly correlated case (U is of the order of the non-interacting bandwidth). Nevertheless, 
our formalism can be used for any other U value belonging to the weak, intermediate, or strong 
coupling regimes. For all methods except the restricted HF (RHF), we have constructed an initial 
guess of the HF transformation such that all symmetries (spin, lattice momentum) are broken. 



This is sometimes referred to as generalized HF (GHF) in the literature 29]. We have converged 
the HF states such that the norm of the gradient is smaller than 10 -5 . For methods involving 
S z projection we have chosen to recover states with S z eigenvalue m = 0, as it is known that at 
half-filling the ground state is always a singlet state (30]. The exact ground state energies, eval- 



2a . have been obtained with an in-house 



uated by solution to the Lieb-Wu equations from Ref. 
Mathematica notebook. 

Table U shows the total energies predicted by a variety of methods for the ground state of the ID 
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TABLE I. Total energies (in units of t, the hopping parameter) for the ground state of the TV-site ID 
Hubbard model Hamiltonian at half-filling with different approximate methods. We have set U — At for all 
calculations. 
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exact£ 
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-1.656854 


-3.748562 


-3.969 123 


A 1 O £2 A rz 
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-4.342 058 


A L* (\ O rr 

-4.6035 
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-2.928 203 


-5.629 064 
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-y.2144 


24 


-6.383016 


-11.258526 


-11.472354 


-11.703 719 


-12.011295 


-13.795 8 


32 


-8.612682 


-15.011368 


-15.224875 


-15.457467 


-15.777256 


-18.3794 


48 


-13.028207 


-22.517052 


-22.730518 


-22.963973 


-23.291156 


-27.5524 


64 


-17.421870 


-30.022 735 


-30.236201 


-30.470037 


-30.798938 


-36.728 7 


96 


-26.187360 


-45.034103 


-45.247569 


-45.481765 


-45.811121 


-55.084 7 


128 


-34.941935 


-60.045471 


-60.258936 


-60.493305 


-60.822554 


-73.4424 


192 


-52.440176 


-90.068206 


-90.281672 


-90.516208 


-90.845 293 


-110.1594 


256 


-69.932961 


-120.090941 


-120.304407 


-120.539025 


-120.868029 


-146.8772 



a Restricted Hartree-Fock, i.e., all symmetries of the Hamiltonian are preserved. 

b Symmetry-broken Hartree-Fock. 

c Non-unitary Hartree-Fock, defined by Eq. 1421 

d S z -projected Hartree-Fock (with S z eigenvalue m = 0). 

e S 2 -projected non-unitary Hartree-Fock, defined by Eq. B3 (with S z eigenvalue m — 0). 
f Obtained by solution to the Lieb-Wu equations of Ref. |28|. 

Hubbard Hamiltonian at half-filling (N = N s , where N is the number of electrons in the system). 
It is evident from the results shown in Table U that nu-HF (defined by Eq. |4"2|) . which uses the 
full flexibility of a non-unitary HF transformation, is able to yield lower energies than standard 
HF. This was expected, since it is at the very least a two-configuration wavefunction. Similarly, 
nu-S 2 HF (defined by Eq. [81~|) yields lower energies than standard S 2 -projected HF. 



31| as the difference with 



It is less evident that the total correlation energy, defined here 
respect to the energy of the broken-symmetry HF solution, should tend to a non-zero constant 
with increasing lattice size. This is the case, as shown in Fig. [TJ In any case, the correlation 
energy per particle predicted by all approximate methods considered in Table U] goes to zero as 
N — > oo. This is a reflection of the limited flexibility that the projected HF and the non- unitary 
based ansatze still have. Note, however, that even if the energy per particle becomes the same as 
N — > oo, the total energy and the wavefunction itself are different from the symmetry broken HF 
solution. 

It is interesting to observe that simple ansatze such as nu-HF or nu-S^HF can be useful to 
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FIG. 1. Total correlation energy (in units of t), predicted by S Z HF and the non-unitary based ansatze, for 
ID Hubbard model calculations as a function of the number of sites N s . The calculations were performed 
at half-filling, with U = At. The correlation energy has been defined with respect to the broken-symmetry 
HF solution. 



describe finite-size lattices where they can capture a significant part of the correlation. For N = 12, 
for which the exact ground state energy is —6. 9204 1, nu-HF recovers 17%, S Z HF recovers 34%, and 
nu-S 2 HF recovers 53 % of the missing correlation energy in the broken-symmetry HF solution. Full 
spin and linear momentum projection may be used to recover even a larger fraction of correlation 
energy, as has been shown for projected HF methods in small size Hubbard ID or 2D lattices 
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Comparison with other two-determinant approaches 



The results shown so far indicate that nu-HF and nu-S 2 HF improve upon HF and S 2 HF, respec- 
tively. This is due to a combination of the more general canonical transformation being used and 
the fact that nu-HF and nu-S 2 HF are explicitly constructed as two-determinant configurations. 

It is interesting to compare the non-unitary based ansatze discussed in this paper with other two- 
determinant ansatze resulting from a single, unitary canonical transformation. We have already 
discussed that more general two-determinant ansatze, where each configuration results from an 
independent HF-transformation, have the same flexibility as the non-unitary approaches considered 
in this work, something we have verified numerically. 

One can think of several ways to construct a two-determinant ansatz based on a single, unitary 
canonical transformation. Our experience shows that symmetry-projection approaches are very 
effective in capturing electron correlations. In this sense, several two-element symmetry groups 
can be used in the ID periodic Hubbard Hamiltonian to build a two-state Goldstone manifold: the 
complex-conjugation group built with the elements {I, K}, where I is the identity operator and K 
is the complex conjugation operator, the time-reversal group built with the elements {I, O}, where 
= exp(i -/r S y ) K is the time-reversal operator, or the C 2 group for even lattices built with the 
elements {/, Cn s / 2 }, where Cjv s /2 is the operator performing a 180-degree rotation of the lattice. 
We here consider the complex conjugation group as a representative example. In this subsection, 
we compare our non- unitary based ansatze with KHF, or complex-conjugation restored HF, and 
KS 2 HF, or complex-conjugation and S 2 -projected HF: 



where \Q) is an iV-particle Slater determinant and P " is the S z projection operator (onto m = 0). 

Figure [2] shows the correlation energy per electron predicted by a variety of approximate methods 
for a 14-site periodic ID Hubbard model as a function of the hole-filling (N/N s ). It is interesting 
to note that at half-filling (N/N s = 1) nu-HF and KHF yield exactly the same correlation energy. 
In this sense, the full flexibility of the non-unitary transformation is not being exploited in the 
solution. At other fillings, on the other hand, nu-HF is able to improve substantially over KHF. In 
contrast, nu-S 2 HF yields lower energies (or larger correlation energies) at all fillings, even though 
the improvement is only marginal in some cases. 

Overall, there is no guarantee that introducing more flexibility into an approximate wavefunction 




(88) 



^KS Z HF 



) = d P s *\<&) + c 2 P Sz k\$) 



(89) 
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FIG. 2. Correlation energy per electron (in units of t), predicted by a variety of approximate methods for 
a 14-site periodic ID Hubbard model as a function of N/N s . The correlation energy has been defined with 
respect to the broken-symmetry HF solution. We were unable to converge KHF for N = 8. 

will result in lower energies for every system. We have shown, however, that ansatze based on a 
non- unitary canonical transformation yield lower energies than HF or projected-HF methods. They 
even yield lower energies than KHF or projected-KHF solutions in some cases, despite the fact that 
complex-conjugation projected wavefunctions are also two-determinant configurations, even if they 
result from a single, unitary canonical transformation. 

IX. CONCLUSIONS 

The HF and the HFB wavefunctions constitute the building blocks upon which more elaborate 
many-body methods rely. They are built out of a set of independent quasi-particles resulting from 
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a linear unitary canonical transformation of elementary fermion operators. In this work, we have 
explored the possibility of relaxing the unitarity condition within a HF-type formalism in order to 
have more variational flexibility in the considered wavef unctions. 

The properties of iV-particle Slater determinants constructed from a set of HF-type operators 
resulting from a non-unitary canonical transformation of fermion operators have been discussed. 
We have derived the corresponding Thouless' theorem for such states, which allowed us to compute 
matrix elements in an efficient way by application of Wick's theorem 

An ansatz based on a single Slater determinant is incapable of utilizing the full flexibility of 
a non-unitary transformation. We have therefore introduced a two-determinant ansatz, defined 
by Eq. H2l where all the degrees of freedom of a HF-type non-unitary transformation are used. 
This, however, is not a limitation of the non-unitary transformation. One could work with other 
more general ansatze used in many-body theory that utilize an iV-particle Slater determinant as a 
starting point. 

Symmetry-breaking is commonly used within a HF formalism to access relevant correlations 
that are otherwise difficult to obtain starting from a symmetry-preserving Slater determinant. In 
this sense, a non-unitary transformation provides additional degrees of freedom that can be used 
in the variational problem. A symmetry-broken wavefunction is, nevertheless, still unphysical; 
we advocate the use of projection techniques out of a symmetry-broken intrinsic state, within a 
variation- after-projection approach, to access the relevant correlations resulting from large quantum 
fluctuations. This can be done, as we have shown in the present work, in combination with a non- 
unitary canonical transformation, affording even more flexibility than that which a projected HF 
state based on a unitary HF transformation has. A non-unitary based projected HF scheme aims to 
provide an accurate description of a many-particle system with a limited number of configurations, 
still a far-reaching problem in fields such as nuclear and condensed matter physics as well as in 
quantum chemistry. 

Finally, we note that our formalism can also be used in the optimization of A r -particle Slater 
determinants that are considered as approximations to the left- and right-eigenvectors of non- 
Hermitian Hamiltonians. In particular, our work can be directly applied to non-Hermitian Hamil- 
tonians with real eigenvalues, such as those resulting from similarity transformations of a standard 
Hermitian one. 

The extension of this work to the full non-unitary Bogoliubov transformation is possible and 
will be presented in a forthcoming publication 
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Appendix A: Proof of Thouless' thorem 

In order to prove the extension to Thouless' theorem stated in section HVl we start by introducing 
the operators {b k , b p , bh, bp} and {d\, d p , dh, dp}, such that {b k , b p } kill the vacuum [<&o) an d \b\, b p } 
kill the vacuum |3>o}> while {d\, d p } annihilate the vacuum and {cft h ,d p } annihilate the vacuum 
|<3?i). We assume that both sets obey the anti-commutation rules defined by Eq. [T4j We explicitly 
write these operators in the form of Eqs. [9] and [T2l that is, 

3 j 

where the superscripts on the matrices indicate the state to which the operators correspond. 

We can now relate the operators {d k ,d p } to the operators {b k ,b p } by using the inverse trans- 
formation discussed in Eq. [15j We arrive at 

4 = E L ^ 6 I' + E y A^ ( Ala ) 

h' V 

d p = Y^ M p , p bp, + Y hp h , (Alb) 

p> h 

where we have set 

L Wh = (D^D 1 ) , (A2a) 

V / h'h 

M plp =(D^D l ) , (A2b) 

v / p'p 

Y ph =(D^D 1 ) , (A2c) 

V / ph 

Y hp =(D^D 1 ) . (A2d) 

V /hp 
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We now assume that the N x N matrix L is invertible, which is only true if ($o|^i) 7^ ( see 

Eq. l2Tjh In such a case, the matrix M is also invertible. We now introduce the operators 

4 = E( L ^)^4- (A3a) 

h' 

d p = Y,[M- X ) plp d pl . (A3b) 

Inserting Eq. IA3I into Eq. IA11 we arrive at 

d[ = b[ + Y, Z Ph~ b l ( A4a ) 
P 

d p = b p + J2 W Phh, (A4b) 



where we have set 



Z P h ~ ^2 Y P*h> ( L * l ) h 'h 
h' 

W ph = J2Yh P ' [M-% p - 



(A5a) 
(A5b) 



In fact, by computing the ant i- commutation rules among the operators {dt, dp}, one can readily 
conclude that W = —Z. This also implies that if L is invertible, then so is M. The transformed 
operators become 



d{ = bl + J2 Z Phbl 
p 

dp = bp — Z p h b~h- 



(A6a) 
(A6b) 



We are now in a position to investigate whether the transformed operators, defined by Eqs. IA6al 
and IA6b"| annihilate the vacuum defined by Eq. [19j We start by evaluating the commutators 



&L ex P ( z~2, Z P' h ' h l' hh ' 

yP'h' 



E Z ph b l eX P E Z P' h ' b l' bh ' 



(A7a) 



p'h' 



bp, ex P ( Z P' h ' b l' bh ' 
p'h' 



^2 Z ph b h J exp ^ Zp> h > bp, by . (A7b) 

\ h J \p'h' I 



The operators from Eqs. IA6al and IA6bl act on the vacuum of Eq. [T9l as 



4 ex P /2 Z P hb l bh l*°> 



Z,h% + Yl Z r>H> ex P \Y. Z W tf/ h' l*o> = 0, (A8f 



ph 



p'h' 



d p exp Z ph b l b h l$o) = Yl Z P hbh - Yl Z P hbh exp \ Y1 Z P' h ' b l' bh ' = °- ( A8b ) 

\ph J \ h h J \p'h' J 
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This essentially completes the proof. {d h ,d p } annihilate the r.h.s. of Eq. IA6bl The operators 
{d h ,d p } that kill the vacuum [<&i) on the l.h.s. of Eq. IA6bl are simple linear combinations of 
{d h , dp}; iV-particle Slater determinants built from either sets of operators are the same up to a 
normalization factor. 



Appendix B: Matrix elements appearing in projected states 

Here, we provide explicit formulas for the matrix elements appearing in the energy expression 
and in the local gradient from the variational ansatz based on projected states. 
The overlap kernels appearing in Eq. [83] are evaluated as 

($|-Re|$) = det N D T R(9) D* , (Bla) 

(l|i?<?|$) = det N D T R{9) D* , (Bib) 

{<S>\R e \$} = det N D T R{9) D* , (Blc) 

(3>\R 9 \$) = det N D T R{9) D* . (Bid) 

The Hamiltonian kernels appearing in Eq. [82] are evaluated in terms of transition density 
matrices as 

($ a \HRg\fy) =Tr ( h p aP f ^ j_ I r aP fm n a0 



Tr ( h P ap (9) + - T ap (9)p ap (9)), (B2) 



T^(9) = Y,(i3\m P ?f(9). (B3) 
ji 

The transition density matrices are in turn given by 

W = { ^h R f } = E D * D ™ + E Aft Z$\e) Dl p , (B4a) 



Pk 



{$\Ro\$) 

I 

($\R e m 



ph 



12 w = m ^ R i m = E Dih °t h + E Dih ^ ( B4b ) 



ph 



= = E Aft + E A, ^V) Dtp, (B4c) 



ph 



&® = m ^ R f } = E Dih D t h + E Aft AV (B4d) 
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Here, 



r("), 
' P h 1 

r(12), 



= E r Ti? w 5 *) pfe , ( £(12) * -1 ^lL' (B5b) 



^£ 2) w = E { DjR ^ D *) h , ( £(22) *" 1 ^)^' ^ B5d ) 



and 



4#(0)=(l3tir(0)2>) A/fc) (B6a) 

•(12), 



^K ) W=(^ t ^W^J^ I (B6b) 

4S(*)=( 5t **W^) fcV (B6c) 
4^(e)=(5t/T(e)5) (B6d) 



The overlap-like matrix elements appearing in the local gradient (Eq. [86]) can be evaluated as 



($\b{b p R e \$) 
(*|^fl|*> 


- V n* 

/ j ^mh 

mn 




Pnm(8) 


($|bj> p ik|$) 


= E ®* mh 

mn 


D n p 


Pnm(8) 




= E ^*mh 
mn 


D n p 


Pnm(0) 


($|&J> p i^|¥) 


= E ^*rah 
mn 


D n p 





(B7a) 
(B7b) 
(B7c) 
(B7d) 
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Similarly, the Hamiltonian-like matrix elements in Eq. [86] can be evaluated as 
mb[b P HR e m (*\HRo\$) 

/ , u mh u np PnrnSP) 



+ E E D *mh D np [h ik + Tj£(e)) pf m {0) (d ni - p%(0)) , (B8a) 



mn ik 



(^\b[b p HR e m 



/ ,, u mh u np Pnm\y) 7=, g [TT 



+ E E ^ A* fa* + r '^)) Pl m (°) i S m - PnM) , (B8b) 



St 



m — / , u mh u np Pnm\° I m 

+ E E ^ fa* + r ^(^)) ^m(^) (<*ni " Pn* (*)) , (B8c) 

mn ifc 

(mem <wi*> 

+ EE^ ^ fa* + T ^)) PfmiP) {S„i - pl\{0)) . (B8d) 
mn ifc 
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